Reservoirs of Stability: Flux Tubes in the Dynamics of Cortical Circuits 
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Triggering a single additional spike in a cerebral cortical neuron was recently demonstrated to cause a cascade 
of extra spikes in the network that is likely to rapidly decorrelate the network's microstate. The mechanisms 
involved in this extreme sensitivity of cortical networks are currently not well understood. Here, we show 
in a minimal model of cortical circuit dynamics that exponential state separation after single spike and even 
single synapse perturbations coexists with dynamical stability to infinitesimal state perturbations. We propose 
a unifying picture of exponentially separating flux tubes enclosing unique stable trajectories composing the 
networks' state spaces. 
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Understanding the dynamical characteristics of cerebral 
cortex networks is fundamental for the understanding of sen- 
sory information processing in the brain. Bottom-up inves- 
tigations of different generic neuronal network models have 
led to a variety of results ranging from stable |fll-l3| to chaotic 
dynamics |4|-|6|]. In top-down attempts to construct classifica- 
tion and discrimination systems with such networks, the 'edge 
of chaos' was proposed to be computationally optimal [7j. 
Near this transition between ordered and chaotic dynamics, 
a network can combine the fading memory and the separa- 
tion property, both of which are important for the efficacy of 
computing applications 1^. While fading memory (informa- 
tion about perturbations of the microstate die out over time) is 
achieved by a stable dynamics, the separation property (dis- 
tinguishable inputs lead to significantly different macrostates) 
is best supported by a chaotic dynamics. 

Widely used in reservoir computing [8] and one of the most 
simple models of cortical circuits are networks of randomly 
coupled inhibitory leaky integrate and fire (LIF) neurons 
These networks exhibit stable chaos, characterized by stable 
dynamics with respect to infinitesimal perturbations despite 
an irregular network activity [Ij 2]. They thus exhibit fading 
memory. Whether and how such networks realize the separa- 
tion property is however unclear 

Motivated by the recent observation that real cortical net- 
works are highly sensitive to single spike perturbations we 
examine in this letter how single spike and single synapse per- 
turbations evolve in a formally stable model of generic corti- 
cal circuits. We show that random networks of inhibitory LIF 
neurons exhibit negative definite Lyapunov spectra, confirm- 
ing the existence of stable chaos. The Lyapunov spectra are 
invariant to the network size, indicating that stable dynamics 
is representative for large networks, extensive and preserved 
in the thermodynamic limit. Remarkably, in the limit of large 
connectivity, perturbations decay as fast as in uncoupled neu- 
rons. Single spike perturbations induce only minute firing rate 
responses but surprisingly lead to exponential state separation 
causing complete decoherence of the networks' microstates 
within milliseconds. By examining the transition from un- 
stable dynamics to stable dynamics for arbitrary perturbation 
size, we derive a picture of tangled flux tubes composing the 



networks' phase space. These flux tubes form reservoirs of 
stability enclosing unique stable trajectories, whereas adja- 
cent trajectories separate exponentially fast. In the thermody- 
namic limit the flux tubes become vanishingly small, implying 
that even in the limit of infinitesimal weak perturbations the 
dynamics would be unstable. This contradicts the prediction 
from the Lyapunov spectrum analysis and reveals that charac- 
terizing the dynamics of such networks qualitatively depends 
on the order in which the weak perturbation limit and the ther- 
modynamic limit are taken. 

We studied large sparse networks of N LIF neurons ar- 
ranged on directed Erdos-Renyi random graphs of mean in- 
degree K. The neurons' membrane potentials Vi G (— oo, Vj) 
with i = 1 . . . N satisfy 



(1) 



between spike events. When Vi reaches the threshold Vt = 
1, neuron i emits a spike and Vi is reset to Vr = 0. The 
membrane time constant is denoted r,„. The synaptic input 
currents are 



hit) = VKh - ^T„ 
V K 



5{t~t^h. (2) 
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composed of constant excitatory external currents v KIq and 
inhibitory nondelayed 5 pulses of strength —Jc/^fK, re- 
ceived at the spike times tj"' of the presynaptic neurons j G 
pre(i). The external currents /q were chosen to obtain a target 
network-averaged firing rate v. 

Equivalent to the voltage representation, Eq. ([T), is a 
phase representation in which each neuron is described by 
a phase 0^ G [— oo, 1], obtained through 0,; = —h\{iVi — 
^/Kh)/iVR - /r/o))/Tf== (where if'- = -ln((V^T - 
VkIq) I (Vr — /q)) is the interspike interval of an isolated 
neuron), with a constant phase velocity and the phase transi- 
tion curve f/((/>0 = - ln(cxp(-0,7f'-) + Jo/(A'/o))/7f"^ 
describing the phase updates at spike reception. The neurons' 
phases thus evolve, from time ts_i after the last spike in the 
network through the next spike time ts, at which neuron j* 
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Figure 1: The balanced state in inhibitoiy LIF networks: (a) Asyn- 
chronous irregular spike pattern of 30 randomly chosen neurons, 
(b) fluctuating voltage trace of one neuron (voltage increased to 
= 2 at spikes), (c),(d) broad distributions of individual neu- 
rons' firing rates v and coefficients of variation cv, (e) network- 
averaged firing rate i> and synchrony measure x versus predicted 
rate i4,ai = Io/{JoT,n) (dotted line: guide to the eye for v = i^bai, 
X = [sTD(^ )] where [•] denotes population average), (parameters: 
iV = 10 000,' A" = 1000, u = 10 Hz, Jo = 1, t„, = 10 ms). 



fires, according to the map 



^free 



(b^its-l) + {ts-ts-l)/Tl' 



fori 7^ i* 
fori = i* 



(3) 

The neurons postsynaptic to the spiking neuron j* are i* E 
post (j*). We used the exact phase map (O for numerically 
exact event-based simulations and to analytically calculate the 
single spike Jacobian Dfi.) = ^J'^*-'^ ■ 



di* (tg) for i = j = i* 

D^J (ts) ={l-d,. for f 
i otherwise. 



(4) 



This matrix depends on the spiking neuron j* and on the 
phases of the spike receiving neurons i* through the derivative 
of the phase transition curve d^. (tg) = d^U{4)i-' (t^)) evalu- 
ated at time just before spike reception |10]. Describing 
the evolution of infinitesimal phase perturbations, the single 
spike Jacobians (|4]i were used for numerically exact calcula- 
tions of all Lyapunov exponents Ai>->Ajvina standard 
reorthogonalization procedure [11]. 

As expected from the construction of the LIF networks, the 
dynamics converged to a balanced state. Figure [1] shows a 
representative spike pattern and voltage trace illustrating the 
irregular and asynchronous firing and strong membrane poten- 
tial fluctuations. A second characteristic feature of balanced 
networks is a substantial heterogeneity in the spike statistics 
across neurons, indicated by broad distributions of coefficients 
of variation (cv) and firing rates ( v). Independent of model de- 
tails, the network- averaged firing rate v in the balanced state 
can be predicted as Pbai ~ lo/i JoTm) flO*]. The good agree- 
ment of this prediction with the numerically obtained firing 
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Figure 2: Stable dynamics with respect to infinitesimal perturba- 
tions: (a) Spectrum of Lyapunov exponents {A^} of networks of 
A'' = 10 000 LIF neurons for different connectivities K, inset: close- 
up of spectra for K = 100 and different network sizes A*', (b),(c) 
largest Lyapunov exponent A2 = Amax and mean Lyapunov exponent 
Amean = X^ili versus Connectivity K and average firing rate 
V (dashed lines: random matrix theory for Amean [TQD, (parameters: 
iV = 100 000, K = 1000, V = 10 Hz, Jo = 1, = 10 ms; 
averages of 10 initial conditions). 



rate confirms the dynamical balance of excitation and inhibi- 
tion in the studied networks. 

Although the voltage trajectory of each neuron and the net- 
work state were very irregular, the collective dynamics of the 
networks was apparently completely stable (Fig. |2]i. For all 
firing rates, coupling strengths and connection probabilities, 
the whole spectrum of Lyapunov exponents (disregarding the 
zero exponent for perturbations tangential to the trajectory) 
was negative, confirming the occurrence of so-called stable 
chaos in LIF networks [1, 2]. The invariance of the Lya- 
punov spectra to the network size TV, to our knowledge, for 
the first time demonstrates that this type of dynamics is ex- 
tensive. With increasing connectivity K all Lyapunov ex- 
ponents approached a constant Ai w — l/r„i. This is de- 
duced from the mean Lyapunov exponent given by Amean ~ 
-l/r,„ + {Vt - {V))/{VKh) + 0[1/K) in random matrix 
approximation and the numerical observation that the largest 
exponent approached Amean in the large K-Wmii (|10] and 
Fig. 12b)). These results suggest that in the thermodynamic 
limit arbitrary weak perturbations decay exponentially on the 
single neuron membrane time constant. As will become clear 
in the following, however, this issue is quite delicate. 

Experimentally realizable and well-controlled state pertur- 
bations to the dynamics of cortical networks are the addition 
or suppression of individual spikes li, 12|. Such minimalis- 
tic neurostimulation can elicit complex behavioral responses 
II12II and can trigger a measurable rate response in intact cor- 
tical networks [6]. We therefore examined how such single 
spike perturbations affected the collective dynamics of our 
networks. Here, the simplest single spike perturbation is the 
suppression of a single spike. Figure |3] illustrates the firing 
rate response if one spike is skipped at f = 0. The miss- 
ing inhibition immediately triggered additional spikes in the 
K postsynaptic neurons such that the network-averaged firing 
rate increased abruptly by 5v ^ Kv/N. Since the induced ex- 
tra spikes inhibited further neurons in the network, the over- 
shoot in the firing rate quickly settled back to the stationary 
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Figure 3: Weak firing rate response after single spike failure: (a) 
Sample spike pattern of 20 randomly chosen neurons (gray: refer- 
ence trajectory, black: single spike skipped at t = 0), (b) network- 
averaged firing rate of reference trajectory v and in response to 
skipped spike v versus time for different connectivities K and net- 
work sizes A'^, (c) number of extra spikes Scxtra = N f{v — i5)dt 
in the entire network versus time (rescaled with average input rate 
Kv), (parameters: v — 10 Hz, Jo — 1, Tm — 10 ms; averages of 
100 initial conditions with 10 000 calculations each). 




Figure 4: Sensitivity to single spike failures: (a) Distance D be- 
tween trajectory after spike failure and reference trajectory versus 
time in log-lin plots for different connectivities K and average fir- 
ing rates (b) pseudo Lyapunov exponent Ap from exponential fits 
D ~ exp(Apt) before reaching saturation versus connectivity K 
and average firing rate v, (c) distance-evolution of all parameter sets 
(rescaled with approximate perturbation strength KJo / VT?) versus 
time (rescaled with average input rate Kv) collapse to characteristic 
exponential state separation with rate Ap ~ 0.9Kv (inset: differ- 
ent network sizes N foi K = 100), (parameters: iV = 100 000, 
K = 1000, V — 10 Hz, Jo — 1, Tm ~ 10 ms; averages of 10 initial 
conditions with 100 calculations each). 



state within a time of order 6t ^ l/{Kv). The overall number 
of additional spikes in the networks therefore was NSvSt w 1 
and the one skipped spike was immediately compensated by a 
single extra spike 1 10]. 

Even though the failure of one individual spike resulted in 
very weak and brief firing rate responses, it nevertheless in- 
duced rapid state decoherence. Figure 2]displays the distance 
D{t) — jj J2i ~ 4>i{t)\ between the perturbed trajec- 

tory (spike failure at < = 0) and the reference trajectory. After 
the spike failure, all trajectories separated exponentially fast at 
a surprisingly high rate. Because this exponential separation 
of nearby trajectories is reminiscent of deterministic chaos, 
we call its separation rate the pseudo Lyapunov exponent Ap. 
The pseudo Lyapunov exponent was network size invariant. 
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Figure 5: (Color online) Sensitivity to finite-size perturbations: (a) 
Distance D between perturbed and reference trajectory measured 
at spike times of reference trajectory (projecting out possible time- 
shifts) for perturbations of strengths e = 0.00002, 0.002, 0.2 in log- 
lin plots (gray lines: 20 examples for initial perturbations of same 
size pointing in different random directions perpendicular to trajec- 
tory, color lines: averages of exponentially separating/converging 
cases), (b) probability Ps of exponential state separation versus per- 
turbation strength e in lin-log plot (dashed line: fit to Psie) = 
1 — exp(— e/eft), dotted line: characteristic perturbation size eft sep- 
arating stable from unstable dynamics, shaded areas: strengths cor- 
responding to single synapse and single spike failures), (c) character- 
istic perturbation size er versus network size A'', connectivity K and 
average firing rate i? in log-log plots, (d) symbolic picture of stable 
flux tubes with radius ett (stable dynamics inside flux tube but expo- 
nential separation of adjacent flux tubes), (parameters: A'^ = 10 000, 
K = 1 000, u — 10 Hz, Jo — 1, Tm = 10 ms; averages of 10 initial 
conditions with 100 calculations and 100 random directions each). 



but showed a completely different behavior compared to the 
classical Lyapunov exponents. With increasing connectivity, 
it appears to diverge linearly Ap ~ KD. It is thus expected to 
grow to infinity in the high connectivity limit, reminiscent of 
binary neuron networks exhibitingan infinite Lyapunov expo- 
nent in the thermodynamic limit ij] . 

In the same balanced LIF networks, we thus find stable dy- 
namics in response to infinitesimal peiturbations and unstable 
dynamics in response to single spike failures. To further ana- 
lyze the transition between these completely opposite behav- 
iors, we applied finite perturbations of variable size perpendic- 
ular to the state trajectory (Fig. |5). Depending on the pertur- 
bation strength e and direction Scj) (with J^i ^'t'l = l)- the per- 
turbed trajectory either converged back to the reference trajec- 
tory or diverged exponentially fast. The probability Pg (e) that 
a perturbation of strength e induced exponential state separa- 
tion was very well-fitted by Ps(£) = 1 — exp(— e/eft). Hence, 
eft is a characteristic phase space distance separating stable 
from unstable dynamics. Intriguingly, this distance decreased 
as Eft ^ 1/ (vKND). For large K and N the dynamics in the 
thermodynamic limit (N — > oo) would be unstable even to 
infinitesimal perturbations (e — > 0). Contrary, the analysis of 
the Lyapunov spectra has shown that taking the limit e ^ 
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first and then N oo yields stable dynamics. Thus, the order 
of the limits appears crucial in defining the dynamical nature 
of balanced LIF networks. 

The evolution of finite perturbations suggests a picture 
of stable flux tubes around unique trajectories (Fig. Eld)). 
Perturbations within these flux tubes decayed exponentially, 
whereas perturbations greater than the typical flux tube ra- 
dius £ft ^ l/{y/KND) induced exponential state separation. 
Single synaptic failures correspond to small perturbations of 
size £syn ~ Jo / vKN and therefore had a N and K inde- 
pendent probability of inducing exponential state separation. 
This probability increased linearly with the average firing rate 

Summarizing, our analysis revealed the co-occurence of 
dynamical stability to infinitesimal state perturbations and 
sensitive dependence on single spike and even single synapse 
perturbations in the dynamics of networks of inhibitory LIF 
neurons in the balanced state. They exhibit a negative defi- 
nite extensive Lyapunov spectrum that at first sight suggests 
a well-defined thermodynamic limit of the network dynamics 
characterized by stable chaos as previously proposed 
In this dynamics, single spike failures induce extremely weak 
firing rate responses that become basically negligible for large 
networks. Nevertheless, such single spike perturbations typi- 
cally put the network state on a very different dynamical path 
that diverges exponentially from the original one. The rate of 
exponential state separation was quantified with the so called 
pseudo Lyapunov exponent Ap. The scaling of Ap ^ KD im- 
plies extremely rapid, practically instantaneous, decorrelation 
of network microstates. Our results suggest that the seemingly 
paradoxical coexistence of local stability and exponential state 
separation reflects the partitioning of the networks' phase 
space into a tangle of flux tubes. States within a flux tube are 
attracted to a unique, dynamically stable trajectory. Different 
flux tubes, however, separate exponentially fast. The decreas- 
ing flux tube radius in the large system limit suggests that an 
unstable dynamics dominates the thermodynamic limit. The 
resulting sensitivity to initial conditions is described by the 
rate of flux tube separation, the pseudo Lyapunov exponent, 
that showed no sign of saturation. These findings suggest that 
the previously reported infinite Lyapunov exponent on the one 
hand iQ] and local stabihty on the other hand ^ |2] resulted 
from the order in which the weak perturbation limit and the 
thermodynamic limit were taken. 

For finite networks, the phase space structure revealed here 
may provide a basis for insensitivity to small perturbations 
(e.g. noise or variations in external inputs) and strong sensitiv- 



ity to larger perturbations. In the context of reservoir comput- 
ing, the flux tube radius defines a border between the fading 
property (variations of initial conditions smaller die out 
exponentially) and the separation property (input variations 
larger Eft cause exponentially separating trajectories). Appli- 
cations of LIF neuron networks in reservoir computing may 
thus strongly benefit if the flux tube structure of the network 
phase space is taken into account. Our results of a very high 
pseudo Lyapunov exponent also reveal that the notion of an 
'edge of chaos' is not applicable in these networks. 
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